To reproduce Panels h and i, you need to download the Processed spatial transcriptomic data from the GEO and then obtain the Starfysh (click here to see how to install) immune deconvolution values by following the python code in the “Multiple-slide immune deconvolution by starfysh” script. Gene signatures to use with Starfysh can be found here: ccRcc_TLS_gene_signatures_starfysh.
#Needed objects from Zenodo:#normal_snRNA.rds #Single nuclei data for the Normal adjacent tissue (NATs)#snRNA.rds #Single nuclei data for the tumoral tissue#Manual_signatures from panel a##Import normal data snlibrary(dplyr)library(ggplot2)library(cowplot)library(harmony)library(Seurat)library(ggplot2)library(cowplot)library(ggsci)library(Rcpp)library(RColorBrewer)library(ggpointdensity)library(symphony)library(viridis)options(stringsAsFactors = F)snRNA_NATs <-readRDS("./normal_snRNA.rds")celltype <-c(brewer.pal(11,"Set3")[c(1,3:8,10:11)],pal_aaas(alpha =0.7)(10)[-2],brewer.pal(8,"Set2"),brewer.pal(11,"Paired"))#[13:1]unique(snRNA_NATs$sample)unique(snRNA_NATs$celltype2)#2 samples##Preprocess the same as the tumoral dataVlnPlot(snRNA_NATs, c("nUMI","nGene","percent.mito"), group.by ="sample",pt.size =0)snRNA_NATs <-SetIdent(snRNA_NATs, cells =NULL, value = snRNA_NATs@meta.data$celltype2)snRNA_NATs <-NormalizeData(snRNA_NATs)snRNA_NATs <-FindVariableFeatures(snRNA_NATs, nfeatures =2000)VariableFeaturePlot(snRNA_NATs)snRNA_NATs <-ScaleData(snRNA_NATs, verbose = T)snRNA_NATs <-RunPCA(snRNA_NATs, npcs =50, verbose =FALSE)PCAPlot(snRNA_NATs,cols=celltype)ElbowPlot(snRNA_NATs, ndims =50)snRNA_NATs <-RunHarmony(snRNA_NATs, c("sample"),max.iter.harmony =20)snRNA_NATs <-FindNeighbors(snRNA_NATs, reduction ="harmony", dims =1:50, do.plot = T)snRNA_NATs <-FindClusters(snRNA_NATs, resolution =1)snRNA_NATs <-RunUMAP(snRNA_NATs, reduction ="harmony", dims =1:50, n.components =2)snRNA_NATs <-CellCycleScoring(snRNA_NATs, g2m.features = cc.genes$g2m.genes,s.features = cc.genes$s.genes)###Take out the cells marked as unknownsnRNA_NATs<-subset(snRNA_NATs,celltype2!="unknown")##Proyect in UMAP expression of GABA + celltypes###Calculate signature for GABA_synthesis etc...geneset<-list("GABA_sinthesis"=Manual_signatures$GABA_synthesis,"GABA_receptor"=Manual_signatures$GABA_receptor)##Hallmark 54 genessnRNA_NATs <-AddModuleScore(object = snRNA_NATs, features =list(Manual_signatures$GABA_synthesis), name ="GABA_synthesis")##GABA_receptorsnRNA_NATs <-AddModuleScore(object = snRNA_NATs, features =list(Manual_signatures$GABA_receptor), name ="GABA_receptor")##DOtplot###Panel b leftDotPlot(snRNA_NATs, features =c("GABA_synthesis1","GABA_receptor1"),cluster.idents = T,group.by ="celltype2") +RotatedAxis() +scale_colour_gradient2(low ="dodgerblue2", mid ="grey90", high ="red",midpoint =0) +ylab("celltype2")+xlab("Signatures")+guides(color =guide_colorbar(title ="Score"))##save RDSsaveRDS(snRNA_NATs,"./snRNA_NATs.rds")###For the tumoral:#combined_matrixsnRNA <-readRDS("./snRNA.rds")##Check object..VlnPlot(snRNA, c("nUMI","nGene","percent.mito"), group.by ="sample",pt.size =0)#snRNA$IMM_subtype<-paste0("IMM_",snRNA$class)###Process datasnRNA <-SetIdent(snRNA, cells =NULL, value = snRNA@meta.data$celltype2)snRNA <-NormalizeData(snRNA)snRNA <-FindVariableFeatures(snRNA, nfeatures =2000)VariableFeaturePlot(snRNA)snRNA <-ScaleData(snRNA, verbose = T)snRNA <-RunPCA(snRNA, npcs =50, verbose =FALSE)celltype <-c(brewer.pal(11,"Set3")[c(1,3:8,10:11)],pal_aaas(alpha =0.7)(10)[-2],brewer.pal(8,"Set2"),brewer.pal(11,"Paired"))#[13:1]PCAPlot(snRNA)ElbowPlot(snRNA, ndims =50)snRNA <-RunHarmony(snRNA, c("sample"),max.iter.harmony =20)snRNA <-FindNeighbors(snRNA, reduction ="harmony", dims =1:50, do.plot = T)snRNA <-FindClusters(snRNA, resolution =1)snRNA <-RunUMAP(snRNA, reduction ="harmony", dims =1:50, n.components =2)snRNA <-CellCycleScoring(snRNA, g2m.features = cc.genes$g2m.genes,s.features = cc.genes$s.genes)###Calculate signaturesDefaultAssay(snRNA) <-"SCT"##Hallmark 54 genessnRNA <-AddModuleScore(object = snRNA, features =list(Manual_signatures$GABA_synthesis), name ="GABA_synthesis")##GABA_receptorsnRNA <-AddModuleScore(object = snRNA, features =list(Manual_signatures$GABA_receptor), name ="GABA_receptor")snRNA$celltype3<-snRNA$celltypesnRNA$celltype3[snRNA$celltype2=="1_B"]<-"Bcells"snRNA$celltype3[snRNA$celltype2=="2_Plasma"]<-"Plasma"##save RDSsaveRDS(snRNA,"./snRNA.rds")
Panel c
Code
#Needed objects from panel b#snRNA processed after running panel blibrary(dplyr)library(ggplot2)library(cowplot)library(harmony)library(Seurat)library(ggplot2)library(cowplot)library(ggsci)library(Rcpp)library(RColorBrewer)library(ggpointdensity)library(symphony)library(viridis)library(clustree)library(cluster)options(stringsAsFactors = F)snRNA <-readRDS("./snRNA.rds")Epi_snRNA <- snRNA[,snRNA$celltype2=="Malignant"]##Epi_snRNA <- readRDS("./Epi_snRNA.rds")##Process dataEpi_snRNA <-NormalizeData(Epi_snRNA)Epi_snRNA <-FindVariableFeatures(Epi_snRNA, nfeatures =2000)VariableFeaturePlot(Epi_snRNA)Epi_snRNA <-ScaleData(Epi_snRNA, verbose = T)#Performed an usupervised clustering changing resolution parametersdr_dim.use<-1:30Epi_snRNA <-RunPCA(object = Epi_snRNA, npcs =30, verbose = T) dims.use = dr_dim.use##Get significant PCAsElbowPlot(Epi_snRNA, ndims =30)print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")#Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.#jackstrawnumber_pcs<-30Epi_snRNA =JackStraw(Epi_snRNA, dims =max(1:number_pcs))Epi_snRNA =ScoreJackStraw(Epi_snRNA, dims =1:number_pcs)dims.use = Epi_snRNA@reductions$pca@jackstraw@overall.p.valuesdims.use = dims.use[dims.use[, 2] <0.05, 1] #Generates a list of the dimensions to use that are the most informative...##Apparently all dimensions 1:30 are significant...Epi_snRNA <-FindNeighbors(Epi_snRNA, dims = dims.use,reduction ="harmony",do.plot = T)#Epi_snRNA <- FindNeighbors(Epi_snRNA, reduction = "harmony", dims = 1:50, do.plot = T)#####Find optimal cluster resulutionprint("Performing Analysis to determine the best resolution parameter for clustering")pcs<-as.matrix(Epi_snRNA@reductions$pca@cell.embeddings)d <-as.dist(1-cor(t(pcs)))res <-seq(0.1, 1.0, by =0.1)res<-c(0.025,0.05,res)library(cluster)clust_nb <-c() v<-c()for (i in1:length(res)) {tryCatch({print(i) Epi_snRNA <-FindClusters(Epi_snRNA, verbose =FALSE, resolution = res[i]) clust_nb <-c(clust_nb, length(levels(Epi_snRNA@meta.data[, c(paste0("RNA_snn_res.", res[i]))])))##Get silhouette value sil <-silhouette(as.numeric(Epi_snRNA$seurat_clusters), d, Fun = mean)gc()v <-c(v, mean(sil[, 3]))}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})}Epi_snRNA@meta.data[] <-lapply(Epi_snRNA@meta.data, function(x) {if (is.factor(x)) as.numeric(as.character(x)) else x})names(clust_nb)<-res names(v)<-res library(clustree)pclust =clustree(Epi_snRNA, prefix ="RNA_snn_res.", node_size =10, node_alpha =0.8)ggsave(pclust , file=file.path("clustree.png"), device="png", width=9, height=12) ###Plot the silhouette results... plot(names(v),v)plot(NULL, NULL, xlim =c(min(as.numeric(clust_nb)), max(as.numeric(clust_nb))), ylim =c(-0.1, 0.5), xlab ="Number of Clusters using malignant Spots", ylab ="")rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col ="#EAE9E9", border =FALSE)grid(col ="white", lty =1, lwd =1.5)points(clust_nb,v, pch =19, col = ggplot2::alpha("black", 0.8), cex =1.5)lines(clust_nb, v,col ="#E51718", lwd =2, lty =4)mtext("Average silhouette width", side =2, line =2, cex =1.5, col ="#224A8D", las =3)axis(side =4, at =seq(0, 1, 0.2), labels =c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0"))invisible(dev.copy2pdf(file ="Silhuette_clusters.pdf", width =5.5, height =5))saveRDS(d,"distance_pca_Epi_snRNA.rds")###According to this the best resolution is: 0.1 which gives rise to 5 clustersrm(d);gc()res_final=0.1Epi_snRNA <-FindClusters(object = Epi_snRNA, resolution = res_final)Epi_snRNA <-RunUMAP(Epi_snRNA, reduction ="harmony", dims = dims.use)##CHeck how it turns using pcaEpi_snRNA <-RunUMAP(Epi_snRNA, reduction ="pca", dims = dims.use,reduction.name ="umap_pca")##Visualizep_tmp<-CellDimPlot(srt = Epi_snRNA, group.by =c("seurat_clusters"),reduction ="UMAP",palcolor =c("#fc772d", "#c6d5a8", "#b48f60", "#a75565", "skyblue")#,theme_use = "theme_blank") ##Note: this UMAP does not have the name of the clusters yet, this was produced after (see panel d)ggsave(filename ="UMAP_clusters_ccRCC.pdf",p_tmp,width =12,height =12)#Save objectsaveRDS(Epi_snRNA,"./Epi_snRNA.rds")
Panel d
Code
#Needed objects from panel c#Epi_snRNA processed after running panel c#Databases_signatures objectlibrary(dplyr)library(ggplot2)library(cowplot)library(harmony)library(Seurat)library(ggplot2)library(cowplot)library(ggsci)library(Rcpp)library(RColorBrewer)library(ggpointdensity)library(symphony)library(viridis)library(clustree)library(cluster)#Epi_snRNA <-readRDS("./Epi_snRNA.rds")###Calculate some needed signatures from Supplementary Table S7Normal_kidney_corpuscle_cells_signatures<-read_excel("Supplementary_tables.xls", sheet ="Table S7",skip =1)Normal_kidney_corpuscle_cells_signatures<-as.list(Normal_kidney_corpuscle_cells_signatures)##Normal_kidney_corpuscle_cells_signatures, Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$TAL), name ="Normal_TAL")##Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$PT), name ="Normal_PT")##Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$DCT), name ="Normal_DCT")###Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$CD), name ="Normal_CD")##Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$CNT), name ="Normal_CNT")###Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Normal_kidney_corpuscle_cells_signatures$DTL), name ="Normal_DTL")###library(SCP)##########DEGs and pathway annotation###Takes around 4 hours....be patientEpi_snRNA <-RunDEtest(srt = Epi_snRNA, group_by ="seurat_clusters", fc.threshold =1, only.pos =FALSE,assay ="RNA",BPPARAM = BiocParallel::MulticoreParam(workers =5))gc()VolcanoPlot(srt = Epi_snRNA, group_by ="seurat_clusters")DEGs_all <- Epi_snRNA@tools$DEtest_seurat_clusters$AllMarkers_wilcox##Export as table DEGs_all<-DEGs_all %>%group_by(group1) %>%arrange(p_val_adj)%>%filter(avg_log2FC>=0) write.table(DEGs_all,file ="DEGs_all_clusters_ccrcc_sn.tab",append = F,sep ="\t",quote = F,row.names = F) DEGs <- DEGs_all[with(DEGs_all, avg_log2FC >=0.2& p_val_adj <0.05), ]table(DEGs$group1)##Annotate pathways..library(dplyr)ht_DE_sn_clusters <-FeatureHeatmap(srt = Epi_snRNA, group.by ="seurat_clusters", features = DEGs$gene, feature_split = DEGs$group1,species ="Homo_sapiens", db =c("GO_BP", "KEGG", "Hallmark"), anno_terms =TRUE, #feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")),height =5, width =4,assay ="RNA")print(ht_DE_sn_clusters$plot) ##Heatmap annotated...###GSEA####Enrichment analysis(GSEA)Epi_snRNA <-RunGSEA(srt = Epi_snRNA, group_by ="seurat_clusters", db =c("GO_BP", "KEGG", "Hallmark"), species ="Homo_sapiens",DE_threshold ="p_val_adj < 0.05",BPPARAM = BiocParallel::MulticoreParam(workers =5))##PlotGSEAPlot(srt = Epi_snRNA, group_by ="seurat_clusters", plot_type ="comparison",db=c("GO_BP", "KEGG","Hallmark"),topTerm =5,padjustCutoff =0.05) ##Cluster 3 is Immune related (probably closer to TLS or immune infiltrated zone...)##Let's check other db##Specific group..group_use = "p_gsea_cluster_0<-GSEAPlot(Epi_snRNA, db =c("GO_BP", "KEGG","Hallmark"), group_by ="seurat_clusters", group_use ="0",topTerm =8) ##For 1p_gsea_cluster_1<-GSEAPlot(Epi_snRNA, db =c("GO_BP", "KEGG","Hallmark"), group_by ="seurat_clusters", group_use ="1",topTerm =8) ##For 2p_gsea_cluster_2<-GSEAPlot(Epi_snRNA, db =c("GO_BP", "KEGG","Hallmark"), group_by ="seurat_clusters", group_use ="2",topTerm =8) ##For 3p_gsea_cluster_3<-GSEAPlot(Epi_snRNA, db =c("GO_BP", "KEGG","Hallmark"), group_by ="seurat_clusters", group_use ="3",topTerm =8) p_gsea_cluster_3##For 4p_gsea_cluster_4<-GSEAPlot(Epi_snRNA, db =c("GO_BP", "KEGG","Hallmark"), group_by ="seurat_clusters", group_use ="4",topTerm =8) p_gsea_cluster_4#O: PT, GABA synthesis, GABA receptor, Omega oxidation, MHCII, Fatty acid oxidation, fatty acid degradation, amino acid metabolism#1: CD, CNT, DCT, Angiogenesis, cell cycle, omega oxidation, cytokine response, TNF, Interferon#2: CNT, little GABA synthesis, cell cycle, FAS_pentose, EMT, Hypoxia, pEMT, p53, Oxidative phosphorylation, EMT#3: TAL, CD, DCT, Motzer_fao_AMPK, this is rather patient-specific#4: GABA receptor, TAL, CD, DCT, Angiogenesis####Final annotation of clusters...Epi_snRNA$annoted_ccrcc_clusters<-ifelse(Epi_snRNA$seurat_clusters=="0","0:PT-like", ifelse(Epi_snRNA$seurat_clusters=="1","1:ProInflammatory",ifelse(Epi_snRNA$seurat_clusters=="2","2:p53-EMT", ifelse(Epi_snRNA$seurat_clusters=="3","3:Patient-specific", "4:PI3K")) ))Epi_snRNA$annoted_ccrcc_clusters<-factor(Epi_snRNA$annoted_ccrcc_clusters,levels =c("0:PT-like","1:ProInflammatory","2:p53-EMT","3:Patient-specific","4:PI3K"))#colors_clusters_ccrcc_sn<-c("#fc772d", "#c6d5a8", "#b48f60", "#a75565", "skyblue")names(colors_clusters_ccrcc_sn)<-levels(Epi_snRNA$annoted_ccrcc_clusters)######Calculate some signatures to plot...####For cluster 0:PT-like#Databases_signaturesDatabases_signatures<-readRDS("Databases_signatures.rds") #Available in the Data folder in this github####For cluster 0:PT-like#GOBP_MONOCARBOXYLIC_ACID_CATABOLIC_PROCESSEpi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$GOBP_MONOCARBOXYLIC_ACID_CATABOLIC_PROCESS), name ="Carboxyic_acid_catabolic")#KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATIONEpi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION), name ="VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION")#KEGG_TRYPTOPHAN_METABOLISMEpi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$KEGG_TRYPTOPHAN_METABOLISM), name ="TRYPTOPHAN_METABOLISM")#HALLMARK_FATTY_ACID_METABOLISMEpi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_FATTY_ACID_METABOLISM), name ="FATTY_ACID_METABOLISM")#####For cluster 1:ProInflammatory Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_TNFA_SIGNALING_VIA_NFKB), name ="TNFA_SIGNALING_VIA_NFKB")#Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_INTERFERON_GAMMA_RESPONSE), name ="IFN_gamma")#Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_IL6_JAK_STAT3_SIGNALING), name ="IL6_JAK_STAT3")#####For cluster 2:p53-EMTEpi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$KEGG_OXIDATIVE_PHOSPHORYLATION), name ="OXIDATIVE_PHOSPHORYLATION")#Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION), name ="Hallmark_EMT")#Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$HALLMARK_P53_PATHWAY), name ="p53_pathway")#####For cluster 3:Patient-specific######For cluster 4:PI3K#Epi_snRNA <-AddModuleScore(object = Epi_snRNA, features =list(Databases_signatures$KEGG_PHOSPHATIDYLINOSITOL_SIGNALING_SYSTEM), name ="PI3K_sginaling")###Final DotPlotsignatures_toPlot_final<-c("GABA_synthesis1","GABA_receptor1","Normal_PT1","Carboxyic_acid_catabolic1","VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION1","TRYPTOPHAN_METABOLISM1","FATTY_ACID_METABOLISM1","TNFA_SIGNALING_VIA_NFKB1","IFN_gamma1","IL6_JAK_STAT31","Normal_CNT1","OXIDATIVE_PHOSPHORYLATION1","Hallmark_EMT1","p53_pathway1","Normal_TAL1","PI3K_sginaling1")dotplot_clusters_sn_ccrcc<-GroupHeatmap(srt = Epi_snRNA,features = signatures_toPlot_final,#group.by =c("annoted_ccrcc_clusters"),#heatmap_palette = "YlOrRd",#cell_annotation = c("RECIST_R_nR"),#cell_annotation_palette = c("Dark2"),show_row_names =TRUE, row_names_side ="left",add_dot =TRUE, add_reticle =TRUE,flip = F, column_title_rot =45, show_column_names =TRUE,row_title_rot =45,column_names_rot =45,height =4.3,width =1.5#cluster_rows = T,#cluster_columns = T# heatmap_palette = "magma"# group_palcolor = list(c("#F98982","#09C1C6")))print(dotplot_clusters_sn_ccrcc$plot)
Panel e and f
Code
##Needed objects from panel d#Epi_snRNA processed after running panel d##Needed objects from panel b#snRNA_NATs processed after running panel blibrary(dplyr)library(ggplot2)library(cowplot)library(harmony)library(Seurat)library(ggplot2)#Epi_snRNA <-readRDS("./Epi_snRNA.rds")snRNA_NATs <-readRDS("./normal_snRNA.rds")##Merge the needed cell togetherPT_normal_ccrcc_integrated<-merge(subset(Epi_snRNA,annoted_ccrcc_clusters!="3:Patient-specific"),subset(snRNA_NATs,celltype2=="PT"))#Integrate using Harmony for faster computing...PT_normal_ccrcc_integrated <-NormalizeData(PT_normal_ccrcc_integrated)PT_normal_ccrcc_integrated <-FindVariableFeatures(PT_normal_ccrcc_integrated, nfeatures =2000)VariableFeaturePlot(PT_normal_ccrcc_integrated)PT_normal_ccrcc_integrated <-ScaleData(PT_normal_ccrcc_integrated, verbose = T)PT_normal_ccrcc_integrated <-RunPCA(PT_normal_ccrcc_integrated, npcs =50, verbose =FALSE)PCAPlot(PT_normal_ccrcc_integrated)ElbowPlot(PT_normal_ccrcc_integrated, ndims =50) PT_normal_ccrcc_integrated <-RunHarmony(PT_normal_ccrcc_integrated, c("sample"),max.iter.harmony =20)PT_normal_ccrcc_integrated <-FindNeighbors(PT_normal_ccrcc_integrated, reduction ="harmony", dims =1:50, do.plot = T)PT_normal_ccrcc_integrated <-FindClusters(PT_normal_ccrcc_integrated, resolution =1)PT_normal_ccrcc_integrated <-RunUMAP(PT_normal_ccrcc_integrated, reduction ="harmony", dims =1:50, n.components =2)##VisualizePT_normal_ccrcc_integrated$annoted_ccrcc_clusters_PT<-ifelse(is.na(PT_normal_ccrcc_integrated$annoted_ccrcc_clusters),"normal-PT",PT_normal_ccrcc_integrated$annoted_ccrcc_clusters)saveRDS(PT_normal_ccrcc_integrated,file ="PT_normal_ccrcc_integrated.RDS")#Remove previous objectsrm(Epi_snRNA,snRNA_NATs)######Perform trayectory analysis using SlingShot from SCPlibrary(SCP)#Trajectory inference###################RunSlingshotPT_normal_ccrcc_integrated <-RunSlingshot(srt = PT_normal_ccrcc_integrated, group.by ="annoted_ccrcc_clusters_PT", reduction ="UMAP",start ="normal-PT")FeatureDimPlot(PT_normal_ccrcc_integrated, features =paste0("Lineage", 1:2), reduction ="UMAP", theme_use ="theme_blank")###THis is the panel eCellDimPlot(PT_normal_ccrcc_integrated, group.by ="annoted_ccrcc_clusters_PT", reduction ="UMAP", lineages =paste0("Lineage", 1:2), lineages_span =0.1)##Dynamic features.... for panel fPT_normal_ccrcc_integrated <-RunDynamicFeatures(srt = PT_normal_ccrcc_integrated, lineages =c("Lineage1", "Lineage2"), n_candidates =200,BPPARAM = BiocParallel::MulticoreParam(workers =15)#,features = genes_to_keep2 #set to genes that are characteristic of each population...along with metaprogrammes (SMP) )DynamicPlot(srt = PT_normal_ccrcc_integrated, lineages =c("Lineage1", "Lineage2"), group.by ="annoted_ccrcc_clusters_PT",features =c("CA9","PAX8","GABA_sinthesis1"),compare_lineages =TRUE, compare_features =FALSE)
Panel g
Code
#You need the spatial transcriptomics and metabolomics from Zenodo:#This is an example for sample R29T#Manual_signatures from Supplementary Table S12spatial_correlation2<-function(spat_obj,var,threshold=0.5,cortest="pearson",colors_pal=c("firebrick1","gray72","palegreen","dodgerblue")){ values <-NULLfor(v in var){if(v %in%rownames(spat_obj)){ spat_obj <-AddMetaData(spat_obj,ifelse(spat_obj@assays$SCT@data[v,] >quantile(spat_obj@assays$SCT@data[v,],threshold),"high","low"),col.name =paste0(v,"_cat")) values[[v]] <- spat_obj@assays$SCT@data[v,] }else{ spat_obj <-AddMetaData(spat_obj,ifelse(spat_obj@meta.data[,v] >quantile(spat_obj@meta.data[,v],threshold),"high","low"),col.name =paste0(v,"_cat")) values[[v]] <- spat_obj@meta.data[,v] } } var_cat <-paste0(var,sep="_cat") spat_obj <-AddMetaData(spat_obj,factor(apply(spat_obj@meta.data[,var_cat],1,function(x) paste0(x,collapse ="_")),levels=rev(c("low_low","low_high","high_low","high_high"))),col.name =paste0(var,collapse ="_")) get_cols <-which(table(spat_obj[[paste0(var,collapse ="_")]][,1])!=0) get_cols <-brewer.pal(4,"Spectral")[get_cols]# p1 <- SpatialDimPlot(spat_obj,# paste0(var,collapse = "_"),# cols = get_cols) p1 <-SpatialPlot(spat_obj,paste0(var,collapse ="_"),image.alpha =0,cols = get_cols) corres <-cor.test(values[[var[1]]], values[[var[2]]],method = cortest,na.action="na.omit") table_cont <-spat_obj@meta.data[,var_cat] order_factor <-rev(c("high_high","high_low","low_high","low_low") ) p2 <-plot_contingency_one_stack(data = table_cont,variables =colnames(table_cont),create_labels = F,order_factor = order_factor,palette =rev(brewer.pal(4,"Spectral"))) p_cor <-ifelse(round(corres$p.value,digits =3) ==0,"p < 10e-16",paste0("p = ",round(corres$p.value,digits =3))) p_chisq <-ifelse(round(p2[[2]]$p.value,digits =3) ==0,"p < 10e-16",paste0("p = ",round(p2[[2]]$p.value,digits =3)))#p1 <- p1 + ggtitle(paste0(var,collapse = " "),# paste0("Correlation = ",round(corres$estimate,digits = 3)," ", p_cor)) p1 <- p1 +ggtitle(paste0(var,collapse =" "),paste0("Chisq.test: ",p_chisq, "\n Rho: ", round(corres$estimate,digits =3)))+scale_fill_manual(values=colors_pal) +DarkTheme() + hide_axis#####Change title to red and bold if significant chisq.test (<0.01)if (p2[[2]]$p.value<=0.01&corres$estimate>0) { p1 <- p1 +theme(plot.subtitle =element_text(color="red", face ="bold")) }#compute correlation plots v1 <-setNames(values[[1]],colnames(spat_obj)) v2 <-setNames(values[[2]],colnames(spat_obj)) corr <-plot_correlation(vect_x =v1,vect_y = v2,labls = var,method = cortest) p3 <- corr[[2]]+theme_classic()return(list("spatial"=p1,"cont_table"=p2,"spatial_cor"=p3)) }library(Cardinal)library(ggpubr)library(ggsci)library(ggrepel)########number II#R29_T, high expressionR29_T_neg <-readImzML(name="R29-T-neg", folder ="spatial_metabolomics/") #Path were you have the downloaded the data###This is for number (II) from panel gimage(R29_T_neg,mz=102.05623, normalize.image="none",colorscale=c(rep("black",2),col.map("jet",n=99)),contrast.enhance="histogram",smooth.image ="gaussian")############Now the remaining numbers using spatial transcriptomicsvisium_directory<-"Path_to_visium"#Path to download data from sample R29_Tseurat_obj <- Seurat::Load10X_Spatial(paste0(visium_directory,"/outs"))#Low-quality / dying cells often exhibit extensive mitochondrial contamination seurat_obj[["percent.mito"]] <-PercentageFeatureSet(seurat_obj, pattern ="^MT.")#VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)# Collect all genes coded on the mitochondrial genome mt.genes <-grep(pattern ="^MT-", x =rownames(seurat_obj), value =TRUE)#remove mt genes genes_to_keep <-setdiff(names(which(Matrix::rowSums(seurat_obj@assays$Spatial$counts )>min.genes)),mt.genes)#Keep genes with count superior to 5#filter cells passing threashold seurat_obj_before<-seurat_obj percent.mito_manual=30 seurat_obj <-subset(seurat_obj,features =genes_to_keep, subset = nFeature_Spatial >300& percent.mito < percent.mito_manual)##Spatial spots featuring more than 30% of mitochondrial genes and less than 300 genes were filtered out, as they identified necrotic or damaged tissue areas which were validated by a pathologist (VV). Genes with counts in less than 5 spatial spots were discarded.cat("Spots removed: ", ncol(seurat_obj_before) -ncol(seurat_obj), "\n")cat("Genes kept: ", length(genes_to_keep),"from",nrow(seurat_obj), "\n") ##Normalize expression data in Seurat..seurat_obj <-NormalizeData(seurat_obj, normalization.method ="LogNormalize", scale.factor =10000)#identify highly variable genesseurat_obj <-FindVariableFeatures(seurat_obj, selection.method ="vst", nfeatures =2000)#Scale data prior to dimension reduction ( stored in pbmc[["RNA"]]@scale.data )all.genes <-rownames(seurat_obj)vars.to.regress="percent.mito"seurat_obj <-ScaleData(seurat_obj, features = all.genes,vars.to.regress = vars.to.regress)###Calculate GABA synthesis and GABA receptor score##Gaba synthesis as meanintersect_tmp<-intersect(Manual_signatures$GABA_synthesis,rownames(seurat_obj))print(length(intersect_tmp))seurat_obj <-AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name ="GABA_synthesis_mean")##Gaba receptor as meanintersect_tmp<-intersect(Manual_signatures$GABA_receptor,rownames(seurat_obj))print(length(intersect_tmp))seurat_obj <-AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name ="GABA_receptor_mean")####Import signatures from Supplementary Table S8 to ccrcc_sn_clusters_signatures##Calculate PT signature from sn using the signature list...intersect_tmp<-intersect(ccrcc_sn_clusters_signatures$`0:PT-like`,rownames(seurat_obj))print(length(intersect_tmp))seurat_obj <-AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name ="ccRCC_PT_like")###################Plot for number (i), H&E imagehe <-SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha =c(0,0)) +NoLegend()###################Plot for number (iii), GABA_synthesis_mean spatially#Use Spatalibrary(SPATA2)###Extra for SPATA2tmp_spata<-SPATA2::asSPATA2(object = seurat_obj,sample_name = tmp_name,image_name ="slice1",spatial_method ="Visium")pl3<-SPATA2::plotSurface(object = tmp_spata,color_by =c("GABA_synthesis_mean"), display_image=T,#alpha_by = "SCGB3A1",smooth = F, normalize =F,#smooth_span = 0.5,pt_size =1.2) + ggplot2::labs(title ="GABA_synthesis_mean")###################Plot for number (iv), ccRCC_PT_like spatiallypl4<-SPATA2::plotSurface(object = tmp_spata,color_by =c("ccRCC_PT_like"), display_image=T,#alpha_by = "SCGB3A1",smooth = F, normalize =F,#smooth_span = 0.5,pt_size =1.2) + ggplot2::labs(title ="ccRCC_PT_like")###############Plot for number (v), the spatial correlationtmp<-spatial_correlation2(spat_obj = seurat_obj,var =c("ccRCC_PT_like","GABA_synthesis_mean"),cortest ="spearman")tmp$spatial
Panel h
Code
###To reproduce this panel you need:#1: Download all the slides visium data from GEO#2: Download the file "Signatures_Starfysh.csv" found in the "Data" folder#3: Obtaine starfysh deconvolution values for each slide by Runnning the python script: "Multiple-slide immune deconvolution by starfysh"#Manual_signatureslibrary(Seurat)process_Visium_seurat<-function(seurat_obj_path,percent.mito_manual=30,min.genes=5,nFeature=300,normalization_method="SCTransform",number_pcs=50,use.jack=T,resolution_clus=1,Calculate_MCPcounter=T,generate_plots=T,plots_path=NULL,marker.test="roc",vars.to.regress="percent.mito",only.pos =TRUE, run_DE=T,...){ raw_data_directory <-seurat_obj_path #This must be the "outs" from SpaceRanger#Load object seurat_obj <- Seurat::Load10X_Spatial(raw_data_directory)#Low-quality / dying cells often exhibit extensive mitochondrial contamination seurat_obj[["percent.mito"]] <-PercentageFeatureSet(seurat_obj, pattern ="^MT.")#VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)# Collect all genes coded on the mitochondrial genome mt.genes <-grep(pattern ="^MT-", x =rownames(seurat_obj), value =TRUE)#remove mt genes genes_to_keep <-setdiff(names(which(Matrix::rowSums(seurat_obj@assays$Spatial$counts )>min.genes)),mt.genes)#Keep genes with count superior to 5#filter cells passing threashold seurat_obj_before<-seurat_obj seurat_obj <-subset(seurat_obj,features =genes_to_keep, subset = nFeature_Spatial > nFeature & percent.mito < percent.mito_manual)##Spatial spots featuring more than 30% of mitochondrial genes and less than 300 genes were filtered out, as they identified necrotic or damaged tissue areas which were validated by a pathologist (VV). Genes with counts in less than 5 spatial spots were discarded.cat("Spots removed: ", ncol(seurat_obj_before) -ncol(seurat_obj), "\n")cat("Genes kept: ", length(genes_to_keep),"from",nrow(seurat_obj), "\n") cat("Performing: ", normalization_method," Normalization","\n") ###Data normalizationif(normalization_method=="SCTransform") {if (use.jack) {cat("You selected SCTransform and JackStraw, the following normalization is only to make JackStraw work\n it is not going to be used for other downstream analyses","\n") #Do this only in case we need the JackStraw since it only supports non-SCT assays seurat_obj_tmp_jack<-seurat_obj seurat_obj_tmp_jack <-NormalizeData(seurat_obj_tmp_jack, normalization.method ="LogNormalize", scale.factor =10000)#identify highly variable genes seurat_obj_tmp_jack <-FindVariableFeatures(seurat_obj_tmp_jack, selection.method ="vst", nfeatures =2000)#Scale data prior to dimension reduction ( stored in pbmc[["RNA"]]@scale.data ) all.genes <-rownames(seurat_obj_tmp_jack) seurat_obj_tmp_jack <-ScaleData(seurat_obj_tmp_jack, features = all.genes,vars.to.regress = vars.to.regress) } seurat_obj <-SCTransform(seurat_obj, assay ="Spatial", verbose = T) }else{#global scaling normalization seurat_obj <-NormalizeData(seurat_obj, normalization.method ="LogNormalize", scale.factor =10000)#identify highly variable genes seurat_obj <-FindVariableFeatures(seurat_obj, selection.method ="vst", nfeatures =2000)#Scale data prior to dimension reduction ( stored in pbmc[["RNA"]]@scale.data ) all.genes <-rownames(seurat_obj) seurat_obj <-ScaleData(seurat_obj, features = all.genes,vars.to.regress = vars.to.regress) }##reduction, calculate PCA, you can manually change the number of pcs to use... seurat_obj <-RunPCA(seurat_obj, features =VariableFeatures(object = seurat_obj),npcs = number_pcs)###Perform further analysis using the best dimensions or manually set dimensionsif (use.jack==TRUE&normalization_method=="SCTransform") {print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")#Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.#jackstraw seurat_obj_tmp_jack <-RunPCA(seurat_obj_tmp_jack, features =VariableFeatures(object = seurat_obj_tmp_jack),npcs = number_pcs) seurat_obj_tmp_jack =JackStraw(seurat_obj_tmp_jack, dims =max(1:number_pcs)) seurat_obj_tmp_jack =ScoreJackStraw(seurat_obj_tmp_jack, dims =1:number_pcs) dims.use = seurat_obj_tmp_jack@reductions$pca@jackstraw@overall.p.values dims.use = dims.use[dims.use[, 2] <0.05, 1] #Generates a list of the dimensions to use that are the most informative... seurat_obj <-FindNeighbors(seurat_obj, dims = dims.use) seurat_obj <-FindClusters(seurat_obj, resolution = resolution_clus)##Calculate also phenograph clusters..#Other clustering method: https://github.com/sararselitsky/FastPG cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,dims.use], num_threads =35) seurat_obj$cluster_phenograph = cluster$communities#Run UMAP seurat_obj <-RunUMAP(seurat_obj, dims = dims.use)###Run TSNE seurat_obj <-RunTSNE(seurat_obj, reduction ="pca", dims = dims.use,do.fast = T, k.seed =10, check_duplicates =FALSE,perplexity =30) }if (use.jack==TRUE&normalization_method!="SCTransform") {print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")#Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.#jackstraw seurat_obj =JackStraw(seurat_obj, dims =max(1:number_pcs)) seurat_obj =ScoreJackStraw(seurat_obj, dims =1:number_pcs) dims.use = seurat_obj@reductions$pca@jackstraw@overall.p.values dims.use = dims.use[dims.use[, 2] <0.05, 1] #Generates a list of the dimensions to use that are the most informative... seurat_obj <-FindNeighbors(seurat_obj, dims = dims.use) seurat_obj <-FindClusters(seurat_obj, resolution = resolution_clus)##Calculate also phenograph clusters..#Other clustering method: https://github.com/sararselitsky/FastPG cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,dims.use], num_threads =35) seurat_obj$cluster_phenograph = cluster$communities#Run UMAP seurat_obj <-RunUMAP(seurat_obj, dims = dims.use)###Run TSNE seurat_obj <-RunTSNE(seurat_obj, reduction ="pca", dims = dims.use,do.fast = T, k.seed =10, check_duplicates =FALSE,perplexity =30) }if (use.jack==FALSE) { seurat_obj <-FindNeighbors(seurat_obj, dims =1:number_pcs) seurat_obj <-FindClusters(seurat_obj, resolution = resolution_clus)##Calculate also phenograph clusters..#Other clustering method: https://github.com/sararselitsky/FastPG cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,1:number_pcs], num_threads =35) seurat_obj$cluster_phenograph = cluster$communities#Run UMAP seurat_obj <-RunUMAP(seurat_obj, dims =1:number_pcs)###Run TSNE seurat_obj <-RunTSNE(seurat_obj, reduction ="pca", dims =1:number_pcs,do.fast = T, k.seed =10, check_duplicates =FALSE,perplexity =30) }if (Calculate_MCPcounter) {#compute mcp scoresif (normalization_method=="SCTransform") { mcp_scores <-MCPcounter.estimate(seurat_obj[["SCT"]]$data,featuresType ="HUGO_symbols") }if (normalization_method!="SCTransform") { mcp_scores <-MCPcounter.estimate(seurat_obj[["Spatial"]]$data,featuresType ="HUGO_symbols") }rownames(mcp_scores) <-make.names(rownames(mcp_scores)) cell_types <-rownames(mcp_scores)for(cell_type in cell_types){ seurat_obj <-AddMetaData(object= seurat_obj,metadata = mcp_scores[cell_type,],col.name = cell_type) } }if (generate_plots) {####Plot The original lame + "nCount_Spatial", "nFeature_Spatial" and the spots removed...##Get the sample name slide<-sub("\\/outs.*", "", raw_data_directory) slide<-sub(".*\\/", "", slide) he <-SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha =c(0,0)) +NoLegend() +ggtitle(slide) to_plot<-c("nCount_Spatial","nFeature_Spatial") removed=list(removed=setdiff(rownames(seurat_obj_before@meta.data),rownames(seurat_obj@meta.data))) spots_removed<-SpatialDimPlot(seurat_obj_before, cells.highlight =removed,facet.highlight = F, ncol =1,pt.size.factor =1.4)+ggtitle(paste0("Spots removed: ", ncol(seurat_obj_before) -ncol(seurat_obj))) pl<-suppressMessages(lapply(to_plot,function(sign){SpatialPlot(seurat_obj,image.alpha =0, pt.size.factor =1.5,features = sign)+scale_fill_viridis_c(option='turbo') +DarkTheme() + hide_axis +theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())+theme(text=element_text(size=10))+theme(text=element_text(face ="bold"))+theme(legend.text=element_text(size=8))+theme(legend.text =element_text(angle =90,margin =margin(8, 0, 0, 0, "pt"))) })) pl <-append(list(he,spots_removed),pl)#save pdf##Get the sample name slide<-sub("\\/outs.*", "", raw_data_directory) slide<-sub(".*\\/", "", slide)pdf(file =paste0(plots_path,"/",slide,"_",Sys.Date(),"_","Quality_spots_features.pdf"), width =20, height =6)print(patchwork::wrap_plots(pl,ncol =4)) dev.off() ###Plot MCP counter populationsif (Calculate_MCPcounter) {##Get the sample name slide<-sub("\\/outs.*", "", raw_data_directory) slide<-sub(".*\\/", "", slide) to_plot <-make.names(c(cell_types)) he <-SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha =c(0,0)) +NoLegend() +ggtitle(slide) unsup_clustering <-suppressMessages(SpatialDimPlot(seurat_obj,label = T,label.size =8)+NoLegend()) pl<-suppressMessages(lapply(to_plot,function(sign){SpatialPlot(seurat_obj,image.alpha =0, pt.size.factor =1.5,features = sign)+scale_fill_viridis_c(option='turbo') +DarkTheme() + hide_axis +theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())+theme(text=element_text(size=10))+theme(text=element_text(face ="bold"))+theme(legend.text=element_text(size=8))#+theme(legend.text = element_text(angle = 90,margin = margin(8, 0, 0, 0, "pt"))) })) pl <-append(list(he,unsup_clustering),pl)pdf(file =paste0(plots_path,"/",slide,"_",Sys.Date(),"_","MCPcounter.pdf"), width =18, height =8)print(patchwork::wrap_plots(pl,ncol =6)) dev.off() } }if (run_DE) {##Find DE markers across clusters seurat_obj_markers <-FindAllMarkers(seurat_obj, min.pct =0.25, logfc.threshold =0.25,test.use = marker.test)return(list(seurat_obj,seurat_obj_markers)) }else{return(seurat_obj) } }library(SCP)#library(SpaCET)library(MCPcounter)#####Importing starfysh deconvolution values and annotating each slide ##############Import metaDatalibrary(readr)Starfysh_integrated_obs_all <-read_csv("Starfysh_integrated_obs_all.csv")##Path were you stored resultsStarfysh_integrated_obs_all$barcode<-Starfysh_integrated_obs_all$...1rownames(Starfysh_integrated_obs_all)<-Starfysh_integrated_obs_all$...1Starfysh_integrated_obs_all<-Starfysh_integrated_obs_all[,-1]Starfysh_integrated_obs_all<-as.data.frame(Starfysh_integrated_obs_all)#########Charging, preprocessing and calculating GABA_synthesis + GABA_receptor for each slideslides_id<-paste0("P_",seq(01,19,by=1)) ##This is the patient id as downloaded from GEOmeta_all<-data.frame()for (i in1:19) {##Create a folder if it does not exist tmp_name<-slides_id[i]print(tmp_name) file.tmp<-paste0("MyoutsPath","/",tmp_name)if(!file.exists(file.tmp)){dir.create(file.path(file.tmp)) }#setwd(file.tmp) setwd(file.tmp) visiumPath<-"Path_to_outsFolder"#Path to visium results Seurat_obj<-process_Visium_seurat(visiumPath,#Path to outs percent.mito_manual=30,min.genes=5,normalization_method="LogNormalize", #Options: SCTransform or NormalizeDatanFeature=nFeature,number_pcs=50,use.jack=T,resolution_clus=0.5,Calculate_MCPcounter=T,generate_plots=F,plots_path="",marker.test="roc",only.pos =TRUE, run_DE=F )##Normalize using SCTransform but set variables to the same number so we can use the AddModuleScore without lossing genes due to variableFeaturesSeurat_obj <-SCTransform(Seurat_obj, verbose = T, #variable.features.n = nfeatures,assay ="Spatial",vars.to.regress="percent.mito")###Add name of sampleSeurat_obj$Sample<-tmp_name##Calculate GABA signaturesintersect_sign <-intersect(Manual_signatures$GABA_synthesis,rownames(Seurat_obj))print(length(intersect_sign))#Seurat_obj <-AddMetaData(Seurat_obj,apply(as.matrix(Seurat_obj[["SCT"]]$data[intersect_sign,]),2,mean),col.name ="GABA_synthesis_mean")##Now for the receptor intersect_sign <-intersect(Manual_signatures$GABA_receptor,rownames(Seurat_obj))print(length(intersect_sign))#Seurat_obj <-AddMetaData(Seurat_obj,apply(as.matrix(Seurat_obj[["SCT"]]$data[intersect_sign,]),2,mean),col.name ="GABA_receptor_mean")####Add starfysh proportions from All_slides_integrated@meta.data###Add other population proportions...Seurat_obj@meta.data$barcode_starfysh<-paste0(gsub("[_].*","",rownames(Seurat_obj@meta.data)),"-",Seurat_obj$Sample)table(Seurat_obj@meta.data$barcode_starfysh%in%Starfysh_integrated_obs_all$barcode)##sigs_starfysh<-colnames(Starfysh_integrated_obs_all)sigs_starfysh<-sigs_starfysh[-1]#Take out barcodefor (j inc(sigs_starfysh)) { Seurat_obj@meta.data[,j]<-Starfysh_integrated_obs_all[match(Seurat_obj$barcode_starfysh,Starfysh_integrated_obs_all$barcode),j]}###Save object...##Save Seurat objectsaveRDS(Seurat_obj,file =paste0("Seurat_obj_",tmp_name,".rds"))###Add metaData only to another objectmeta_tmp<-Seurat_obj@meta.datameta_all<-rbind(meta_all,meta_tmp)rm(Seurat_obj)gc() }#Save meta_allsaveRDS(meta_all,file ="meta_all.rds")########Calculations for panel hnnls1=c("ccRCC_PT_like_starfysh","ccRCC_Proflammatory_like_starfysh","ccRCC_p53_EMT_starfysh","ccRCC_PI3K_starfysh",sigs_starfysh[1:20])ma=c("GABA_synthesis_mean","GABA_receptor_mean")##Correlation signatures...meta<-meta_allval_all =Reduce(rbind,lapply(ma, function(m){ #ma for modulesReduce(rbind,lapply(c(nnls1), function(ct){ ##Populations res =sapply(unique(meta$Sample), function(o){ val =NAtryCatch(expr = { ct_nei =paste0(ct) ##Recover distance sub = meta[meta$Sample == o ,] ##subset to sample model =lm(sub[,m] ~ sub[,ct_nei]) #module ~ pop_dist est =summary(model)$coefficients[2,1] pval =-log10(summary(model)$coefficients[2,4]) #pvalue val = pval*sign(est)return(val) }, error =function(e){return(NA)}) })#res[abs(res)<pval_thresh] = 0#res = sign(res)#med = sum(res, na.rm = TRUE)#amp = sum(res == sign(med), na.rm = TRUE) med =median(res, na.rm =TRUE) amp =mean(sign(res) ==sign(med), na.rm =TRUE)return(data.frame('module'= m, 'ct'= ct, 'med'= med, 'amp'= amp))return(data.frame('module'= m, 'ct'= ct, 'med'= med, 'amp'= amp)) }))}))head(val_all)###Med is the median of the (+-log10(pvalue)); amp(fraction samples)##Negavite values for med means the correlation is negative;#val_R$module = factor(val_R$module, levels = rev(ma[c(5:8,18,19,15,10,17,16,1:3,9,4,11:14,20)])) ##ordenate according to desire value...val_all$ct =factor(val_all$ct, levels =rev(nnls1)) val_all$module=factor(val_all$module,levels =c("GABA_synthesis_mean","GABA_receptor_mean"))val_all2<-val_allval_all$highlight = (val_all$amp >=0.5) & (abs(val_all$med) >1.3)#Equivalent for pval<0.05val_all$med[val_all$med >4] =4val_all$med[val_all$med <-4] =-4p2R =ggplot(val_all, aes(x=ct,y=module)) +geom_point(shape=21,aes(size=amp,fill=med,colour=highlight)) +scale_fill_gradient2(low =brewer_pal(palette ='RdBu', direction =-1)(n =3)[1],mid =brewer_pal(palette ='RdBu', direction =-1)(n =3)[2], high =brewer_pal(palette ='RdBu', direction =-1)(n =3)[3]) +scale_colour_manual(breaks=c('FALSE','TRUE'),values=c('white','black')) +theme_bw() +theme(panel.grid.minor =element_blank(),panel.grid.major =element_blank(),axis.text =element_text(size=14, colour ="black"),axis.text.y =element_text(size=12, colour ="black"),axis.text.x =element_text(size=12, colour ="black", angle =90, hjust =1),axis.title=element_blank(),panel.border =element_rect(size =0.7, linetype ="solid", colour ="black"))print(p2R+coord_flip())p2R<-p2R+coord_flip()###############Plot boxplot distribution of pvalues ###########For synthesis m<-ma[1] val =Reduce(rbind,lapply(nnls1, function(ct){ res =sapply(unique(meta$Sample), function(o){ val =NAtryCatch(expr = {# ct_dist = paste0(ct, '_dist') sub = meta[meta$Sample == o,] model =lm(sub[,m] ~ sub[,ct]) est =summary(model)$coefficients[2,1] pval =-log10(summary(model)$coefficients[2,4]) pval[pval >10] =10 val = pval*sign(est)return(val) }, error =function(e){return(NA)}) })return(data.frame('Sample'=unique(meta$Sample), 'ct'= ct, 'res'= res)) })) val$ct =factor(val$ct, levels =rev(nnls1)) val$Sample =factor(val$Sample)#val$RECIST_R_nR<-meta$RECIST_R_nR[match(val$Sample,meta$Sample)] p =ggplot(val, aes(x=ct,y=res)) +geom_boxplot(coef =0, outlier.shape =NA) +#geom_violin() + geom_jitter(colour ="gray") +geom_hline(yintercept =log10(0.05), linetype =2) +geom_hline(yintercept =-log10(0.05), linetype =2) +#scale_colour_manual(values = colors, breaks = names(colors)) +ggtitle(m) +theme_bw() +theme(panel.grid.minor =element_blank(),panel.grid.major =element_blank(),axis.text =element_text(size=14, colour ="black"),axis.text.y =element_text(size=12, colour ="black"),axis.text.x =element_text(size=12, colour ="black", angle =90, hjust =1),axis.title=element_blank()) +NoLegend()print(p+coord_flip())p_synthesis<-p+coord_flip()#For receptor m<-ma[2] val =Reduce(rbind,lapply(nnls1, function(ct){ res =sapply(unique(meta$Sample), function(o){ val =NAtryCatch(expr = {# ct_dist = paste0(ct, '_dist') sub = meta[meta$Sample == o,] model =lm(sub[,m] ~ sub[,ct]) est =summary(model)$coefficients[2,1] pval =-log10(summary(model)$coefficients[2,4]) pval[pval >10] =10 val = pval*sign(est)return(val) }, error =function(e){return(NA)}) })return(data.frame('Sample'=unique(meta$Sample), 'ct'= ct, 'res'= res)) })) val$ct =factor(val$ct, levels =rev(nnls1)) val$Sample =factor(val$Sample)#val$RECIST_R_nR<-meta$RECIST_R_nR[match(val$Sample,meta$Sample)] p =ggplot(val, aes(x=ct,y=res)) +geom_boxplot(coef =0, outlier.shape =NA) +#geom_violin() + geom_jitter(colour ="gray") +geom_hline(yintercept =log10(0.05), linetype =2) +geom_hline(yintercept =-log10(0.05), linetype =2) +#scale_colour_manual(values = colors, breaks = names(colors)) +ggtitle(m) +theme_bw() +theme(panel.grid.minor =element_blank(),panel.grid.major =element_blank(),axis.text =element_text(size=14, colour ="black"),axis.text.y =element_text(size=12, colour ="black"),axis.text.x =element_text(size=12, colour ="black", angle =90, hjust =1),axis.title=element_blank()) +NoLegend()print(p+coord_flip())p_receptor<-p+coord_flip()+scale_x_discrete(position ="top")#Plot together...########################Plotting panel hggsave(ggarrange(p_synthesis+ggtitle(NULL),p2R+theme(axis.text.x =element_blank()),p_receptor+ggtitle(NULL),nrow =1,align ="hv"),filename ="Figure3h.pdf",width =14,height =8)
Panel i
Code
##You need to run panel h to reproduce this figure###This is an example for P_6i=6library(Seurat)library(SPATA2)spatial_correlation2<-function(spat_obj,var,threshold=0.5,cortest="pearson",colors_pal=c("firebrick1","gray72","palegreen","dodgerblue")){ values <-NULLfor(v in var){if(v %in%rownames(spat_obj)){ spat_obj <-AddMetaData(spat_obj,ifelse(spat_obj@assays$SCT@data[v,] >quantile(spat_obj@assays$SCT@data[v,],threshold),"high","low"),col.name =paste0(v,"_cat")) values[[v]] <- spat_obj@assays$SCT@data[v,] }else{ spat_obj <-AddMetaData(spat_obj,ifelse(spat_obj@meta.data[,v] >quantile(spat_obj@meta.data[,v],threshold),"high","low"),col.name =paste0(v,"_cat")) values[[v]] <- spat_obj@meta.data[,v] } } var_cat <-paste0(var,sep="_cat") spat_obj <-AddMetaData(spat_obj,factor(apply(spat_obj@meta.data[,var_cat],1,function(x) paste0(x,collapse ="_")),levels=rev(c("low_low","low_high","high_low","high_high"))),col.name =paste0(var,collapse ="_")) get_cols <-which(table(spat_obj[[paste0(var,collapse ="_")]][,1])!=0) get_cols <-brewer.pal(4,"Spectral")[get_cols]# p1 <- SpatialDimPlot(spat_obj,# paste0(var,collapse = "_"),# cols = get_cols) p1 <-SpatialPlot(spat_obj,paste0(var,collapse ="_"),image.alpha =0,cols = get_cols) corres <-cor.test(values[[var[1]]], values[[var[2]]],method = cortest,na.action="na.omit") table_cont <-spat_obj@meta.data[,var_cat] order_factor <-rev(c("high_high","high_low","low_high","low_low") ) p2 <-plot_contingency_one_stack(data = table_cont,variables =colnames(table_cont),create_labels = F,order_factor = order_factor,palette =rev(brewer.pal(4,"Spectral"))) p_cor <-ifelse(round(corres$p.value,digits =3) ==0,"p < 10e-16",paste0("p = ",round(corres$p.value,digits =3))) p_chisq <-ifelse(round(p2[[2]]$p.value,digits =3) ==0,"p < 10e-16",paste0("p = ",round(p2[[2]]$p.value,digits =3)))#p1 <- p1 + ggtitle(paste0(var,collapse = " "),# paste0("Correlation = ",round(corres$estimate,digits = 3)," ", p_cor)) p1 <- p1 +ggtitle(paste0(var,collapse =" "),paste0("Chisq.test: ",p_chisq, "\n Rho: ", round(corres$estimate,digits =3)))+scale_fill_manual(values=colors_pal) +DarkTheme() + hide_axis#####Change title to red and bold if significant chisq.test (<0.01)if (p2[[2]]$p.value<=0.01&corres$estimate>0) { p1 <- p1 +theme(plot.subtitle =element_text(color="red", face ="bold")) }#compute correlation plots v1 <-setNames(values[[1]],colnames(spat_obj)) v2 <-setNames(values[[2]],colnames(spat_obj)) corr <-plot_correlation(vect_x =v1,vect_y = v2,labls = var,method = cortest) p3 <- corr[[2]]+theme_classic()return(list("spatial"=p1,"cont_table"=p2,"spatial_cor"=p3)) } Seurat_obj<-readRDS(file =paste0("Seurat_obj_",tmp_name,".rds"))#For number (i) he_p<-SpatialPlot(Seurat_obj,repel = F,label = F,image.alpha=1,alpha =c(0,0)) +NoLegend() +ggtitle(tmp_name)####Signatures###Extra for SPATA2tmp_spata<-SPATA2::asSPATA2(object = Seurat_obj,sample_name = tmp_name,image_name ="slice1",spatial_method ="Visium")##For number (iii, iv, v, vi, vii, viii, ix)#Returnpl2<-suppressMessages(lapply(c("GABA_synthesis_mean","ccRCC_PT_like_starfysh","ccRCC_Proflammatory_like_starfysh","ccRCC_p53_EMT_starfysh","ccRCC_PI3K_starfysh","GABA_receptor_mean","Starfysh_Macrophage M2"),function(sign){ SPATA2::plotSurface(object = tmp_spata,color_by =c(sign), display_image=T,#alpha_by = "SCGB3A1",smooth = T, normalize =F,smooth_span =0.6,pt_size =1.2) + ggplot2::labs(title = sign) }))#####for numbers (x, xi)##Correlations...var=list(c("GABA_synthesis_mean","ccRCC_PT_like_starfysh"),c("GABA_receptor_mean","Starfysh_Macrophage M2"))##Grab Legend...to add at the end of each graphtmp<-spatial_correlation2(spat_obj = Seurat_obj,var =c("GABA_synthesis_mean","ccRCC_PT_like_starfysh"),cortest ="spearman")tmp$spatial ##Extract legendlibrary(ggpubr)legend_spatial<-get_legend(tmp$spatial+theme(legend.title =element_blank()) )tmpall_plots_cors<-lapply(var,function(v){tryCatch({ res <-spatial_correlation2(spat_obj = Seurat_obj,var = v,cortest ="spearman") res$spatial <- res$spatial + hide_axisif(identical(v,var[length(var)][[1]])){ res$spatial+NoLegend() }else{ res$spatial+NoLegend() } }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})})names(all_plots_cors)<-lapply(var, function(x){paste0(x[1],x[2])})all_plots_cors